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We have obtained by Monte Carlo NVT simulations the constant-volume excess heat capacity of 
square-well fluids for several temperatures, densities and potential widths. Heat capacity is a ther- 
modynamic property much more sensitive to the accuracy of a theory than other thermodynamic 
quantities, such as the compressibility factor. This is illustrated by comparing the reported simu- 
lation data for the heat capacity with the theoretical predictions given by the Barker-Henderson 
perturbation theory as well as with those given by a non-perturbative theoretical model based on 
Baxter's solution of the Percus-Yevick integral equation for sticky hard spheres. Both theories give 
accurate predictions for the equation of state. By contrast, it is found that the Barker-Henderson 
theory strongly underestimates the excess heat capacity for low to moderate temperatures, whereas 
a much better agreement between theory and simulation is achieved with the non-perturbative the- 
oretical model, particularly for small well widths, although the accuracy of the latter worsens for 
high densities and low temperatures, as the well width increases. 



I. INTRODUCTION 

Thermodynamic and structural properties of square-well (SW) fluids have been profusely studied both from theory 
and from computer simulation p|-|38j. From the theoretical side, the first few virial coefficients have been obtained 
QiH>ll3 ^^'^ the radial distribution function has been evaluated from numerical solutions of integral equation theories, 
such as Percus-Yevick H H IM Hp , Yvon-Born-Green [Tl, HNC [l^, MSA [llj33, Rogers-Young ORPA 
[28l |. and HRT [s^- Simpler analytical approximations have also been proposed Il2. 19, ,2jX 22e ^M- The 
thermodynamic properties have been derived from the theoretical structure functions as well as from perturbation 
theory [jj 0, 0, 0j HE 113 ■ Access to the "experimental" properties of SW fluids has been made possible via 
molecular dynamics and Monte Carlo simulations Jj^ JT ^ 1 ^ ll^ IT7 [ flS.. 26 , ..36, 38 1. Special attention has received 



> ■4D. ..J p, oo|. Special attentic 
the determination of the critical point of SW fluids |d, l2l IE ll^ HE illM S S H El El , both from the 
theoretical and simulational viewpoints. The main reason for this wide interest lies in the fact that a SW fluid is 
perhaps the simplest one whose particles have attractive as well as repulsive interactions. In general, theories are 
easier to apply to SW fluids than to other fluids with more realistic potentials. In addition, the SW potential seems 
to be particularly sensitive to the performance of a theory. Therefore, this kind of fluid is an excellent testing-ground 
for many theories of fluids and so the study of SW fluids can be considered as a flrst step towards our understanding 
of the properties of fluids with more sophisticated interactions. There is an additional reason explaining the recent 
revival of interest in SW fluids. The SW potential possesses, besides the diameter of the hard core and the depth of 
the well, an additional parameter measuring the width of the well. This makes the SW pot ential with a small width 
especially suited to model the effective interactions among colloidal particles I2M l3(A ISll l3q . In this context, the 
glass transition and a solid-to-solid isostructural transition 24] have been studied for narrow SW systems. 

Despite the extensive number of studies devoted to the SW fluid, relatively little attention has been paid to several 
thermodynamic properties. This is the case for the heat capacity. To the best of our knowledge, there are available 
d only a few simulation data of this property for SW fluids. Theoretical calculations of the same quantity are 
equally scarce |lC| . In the present paper we have carried out Monte Carlo simulations of the constant-volume excess 
heat capacity Cy of SW fluids for several values of the potential width and, for each of them, for several densities 
and temperatures. Moreover, in order to put clear the sensitivity of this property to the accuracy of a theory, the 
simulation data are compared with the results obtained from the Barker-Henderson (BH) 23| perturbation theory 
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and with those derived from the theoretical model proposed by Yuste and Santos [2j|, recently simplified by Acedo 
and Santos |3^ . 

The paper is organized as follows. In the next section, we summarize the theoretical foundations of the MC procedure 
used and we describe the simulations performed and the results obtained. In Section [llll we present an outline of the 
above-mentioned theories. Finally, in the last section the theoretical results are compared with simulation data and 
discussed. 



II. MONTE CARLO SIMULATIONS 



In a square-well (SW) fluid, particles interact by means of a potential of the form 

{oo if r < a 
~e if a < r < Xa (1) 
if r>Xa 

where A is the potential width in units of the particle diameter a and e is the potential depth. 

Constant- volume averaged excess heat capacity per particle in a SW fluid can be expressed in the form 



C0 _ 1 ((M-(M)f) 
Nk ~ r*2 N 



(2) 



where A'^ is the number of particles, k is the Boltzmann constant, T* — kT/e is the reduced temperature, and M is 
the number of pairs of interacting particles, that is, the number of pairs of particles whose centers lie separated by a 
reduced distance r* = r /a < X. 

The averages involved in equation Q can be calculated by Monte Carlo (MC) simulations in the NVT ensemble. 
Therefore, we have proceeded to calculate by means of MC NVT the constant-volume averaged excess heat capacity 
per particle for SW fluids with well widths A ranging from 1.1 to 1.5. For each value of A, Cy has been evaluated for 
several densities along isotherms. To this end, a system consisting of 512 particles placed in a cubic box with periodic 
boundary conditions was used. Particles were initially placed in a regular configuration and then the system was 
allowed to equilibrate for 2 x 10^ cycles, each of them consisting of an attempt move per particle, the first 10"* cycles 
at a very high temperature and the remaining ones at the desired temperature. The calculation of Cy was performed 
by averaging over the next 5 x 10^ cycles, performing partial averages every 10^ cycles with the aim of estimating the 
statistical error from the standard deviation. The use of such a huge number of cycles in the calculations was motivated 
by the need of ensuring that the values of Cy converged to a constant value, apart from statistical fluctuations. In 
fact, we realized that for low values of the number of cycles used in the calculations, the values of Cy increase with 
the number of cycles used. 

The results are shown in Table|lJ We have considered four isotherms for A — 1.1, 1.2, 1.3 and three isotherms for A — 
1.5. The lowest temperature in each case is larger than the estimated critical temperature [T3Lin . l23 . l25ll2^|33Ll33 . l35| : 
T* ~ 0.5, 0.6, 0.8, 1.2 for A = 1.1, 1.2, 1.3, 1.5, respectively. 



III. THEORY 



A. Barker— Henderson perturbation theory 



In the second-order Barker-Henderson perturbation theory |39l |40j , the free energy is expressed in the form 

F Fo , Fi 1 F2 1 , , 



NkT NkT NkTT* NkT T* 

where Fq is the free energy of the hard-sphere (HS) reference system and Fi and F2 are the first- and second-order 
perturbative terms, respectively. According to this theory, the constant-volume excess heat capacity per particle is 
given by 

^ = (4) 
Nk T*^ NkT' ^ ' 
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where 

^^^^{§p) / [u*i{r)f 9o{r)r^dr (5) 



NkT y^. 

in the so-called macroscopic compressibility approximation, whereas 



= -TrpfcT / [ui (r)] <^ — \ r dr (6) 



iVfcT Jo I 9P 







in the so-called local compressibility approximation. In Eqs. (|SJ| and ®i P = N/V is the number density, = 
ui{r)/e is the perturbative contribution to the potential function, which in a SW potential is u*(r) = — 1 for a < r < 
Act, P is the pressure, and go{r) is the radial distribution function (r.d.f.) of the hard-sphere reference fluid. 

In recent years, several analytical and very accurate expressions for the r.d.f. go{r) of the HS fluid have been 
proposed |4l|, |42, 113 ■ They can be used to determine F2 in expressions © and 10. Regarding {dp/dP)o, which 
appears explicitly in expression Q and implicitly in jnj, it can be obtained from the well-known Carnahan-Starling 
[4JI equation of state 

PoV _ 1 + 7? + 77^ - 7?3 

^° " Mr " {{^^ ' 

where 77 — {7r/6)pa^ is the packing fraction. 

B. Yuste—Acedo— Santos model 



The internal energy can be obtained from the r.d.f g{r) through the energy equation 



U = ^NkT + 2-KNp / u{r)g (r) r^dr, (8) 







whence 

In the special case of the SW potential (^1, Eqs. ^ and © become 

U ^-NkT -UNer] [ g (r*) r*'^dr* , (10) 
2 Ji 





\dg{r*^ 







r*^dr*, (11) 



respectively. The r.d.f g{r*) of the SW fluid depends on the packing fraction 77, the reduced temperature T* and, 
parametrically, on the well width A. In principle, one has to resort to numerical solutions of integral equation theories. 
On the other hand, particularly suitable for the purpose of obtaining the heat capacity is the heuristic model proposed 
by Yuste and Santos [23| and recently simplified by Acedo and Santos which is analytical and fairly accurate. 
Henceforth we will refer to this model as the Yuste- Acedo-Santos (YAS) model. It is based on expressing the Laplace 
transform G{t) of r*g{r*) in the form 



F{t)e ^ 

n=l 

where the auxiliary function F{t) is assumed to have the form [22l l3 



G W = \^u%)e-^ - (-12.)"-^ t [F itr e-, (12) 



1 e'/^' + K,t - (e^/^' - I + K2t) e-(^^~^^' 
1 + Sit + S2t^ + S3t^ ■ 
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The coefficients ifi, K2, Si, S2 and S3 are explicit functions of 77, T* and A determined from consistency conditions. 
We refer the interested reader to Refs. [12,1^1 for further details. The YAS model l(T^ reduces to the exact solutions 
of the Percus-Yevick (PY) equation in the limit of hard spheres (A — > 1 or T* ^ 00) 0,113 j as well as in the limit 
of sticky hard spheres (A ^ 1 and T* ^ with T* ^ — 1/ ln(A — 1)) H?]. From that point of view, the approximation 
(|13|l can be seen as a simple extension to finite widths of Baxter's solution of the PY equation for sticky hard spheres. 
Upon Laplace inversion of Eq. H12(l . the final expression of the r.d.f. reads 

oc 

9 in = r*-' J2 (-12??)""' fn {r* - n) 9 (r* - n), (14) 

n=l 

where the functions fn{f*) are the inverse Laplace transforms of and 8(r*) is Heaviside's step function. 

Therefore, to determine the r.d.f. for r* < n + 1 only the first n terms in the summation H14(l are needed. In 
particular, for the values of A < 2 considered in this paper, one has 



_i 3 



I— I * 

where Zi [i = 1, 2, 3) arc the three roots of the cubic equation 1 + Sit + S2t^ + S^t^ = 0. Inserting Eq. (fT3|l into Eq. 
((TT|l . we finally get 



Nk dT* ^ Si+ 2S2Z^ + 3S3Z^ 



zri_i + (A-zri)e-(^-^) . (16) 



The heat capacity can also be obtained from the YAS r.d.f. by following the virial and compressibility routes to 
the equation of state. The reason for the choice of the energy route ^ is two-fold. First, it is obviously the most 
direct route to determine the heat capacity. Second, we have checked that the other routes yield results that present 
larger deviations from the simulation data. This latter observation is consistent with the case of the PY theory for 
sticky hard spheres ^ and for SW fluids 0, . 



IV. RESULTS AND DISCUSSION 



Results obtained for Cy from the second-order Barker-Henderson perturbation theory within the local compress- 
ibility approximation as well as within the macroscopic compressibility approximation are compared in Figs. with 
the simulation data of Table ^ We can see that although the local compressibility approximation provides a better 
agreement with simulation data, both approximations are rather poor at low temperatures. This might be due either 
to the insufficient accuracy of both the local compressibility and the macroscopic compressibility approximations or 
to the fact that higher order terms, beyond the second one, in the expansion of the Helmholtz free energy in power 
series of the inverse of the reduced temperature, have a nonnegligible contribution to the heat capacity. In order to 
determine which of these two possibilities is the right one, we can use for F2 in Eq. ^ simulation data, thus avoiding 
theoretical approximations. These simulations were performed by Barker and Henderson jl^ who reported the results 
in terms of a function depending on forty five parameters for each density. These parameters were determined from 
a least squares fitting of their simulation data. Since the use of that fitting is somewhat tedious, we have preferred 
to use directly simulation data for F2, which are available for several densities and well widths [53, to determine Cy 
from Eq. As one can see in Figs. results thus obtained are much closer to the theoretical results derived 

from the BH second-order perturbation theory than to the values of Cy obtained from direct simulations, except in 
the high density limit. This means that the main reason of the failure of the Barker-Henderson perturbation theory 
in predicting the heat capacity of SW fluids arises in the truncation of the perturbative series at the level of the 
second order term, the higher order terms having a nonnegligible contribution. This is in contrast with the situation 
for the equation of state [s^ Isij . which is accurately given by the second-order BH perturbation theory even at 
relatively low temperatures for wide ranges of densities and potential wells. The reason is that, as pointed before, 
the constant-volume excess heat capacity is a thermodynamic property particularly sensitive to the performance of a 
theory and therefore the influence of higher order terms, which is small in the equation of state, may be important in 
the heat capacity. This is particularly true for low values of the potential width, since the lower the potential width, 
the slower the convergence of the BH perturbation theory at low temperatures |15) . 

A much better agreement is obtained with YAS theory, Eq. (|16|l . at low to moderate densities, as shown in the 
same figures. This theory, in contrast to the BH theory, provides a better agreement with the simulation data of Cy 
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as the potential width decreases. This is consistent with the fact that, as said before, the YAS model is an extension 
to A > 1 of the PY solution for sticky hard spheres and hence it is expected to be as accurate as the PY theory at 
least for small A — 1. The structural properties predicted by the YAS model for the SW fluid also exhibits a good 
agreement with simulation data for low values of A — 1 whereas the accuracy worsens as A increases |^ . Figures 
n^show that, given a well width A, the YAS values of Cy are more accurate as the temperature increases and/or 
the density decreases. 

In conclusion, we have performed Monte Carlo simulations of the constant-volume excess heat capacity of SW 
fluids of variable width for a wide range of densities and at several characteristic temperatures. This thermodynamic 
quantity vanishes for hard spheres and so it represents an important measure of the influence of attractive forces on 
the state of the fluid. Moreover, the heat capacity seems to provide a rather stringent test to assess the accuracy of 
theoretical app roaches. In this paper we have compared the simulation data with the Barker-Henderson perturbation 
theory [3^l4Clj| and with a non-perturbative theory developed by Yuste, Acedo, and Santos "25, 's?]. While the former 
theory presents a poor performance, which can be attributed to the truncation of the perturvative series to second 
order rather than the inaccuracy of the theory itself, the non-perturbative theory does a fairly good job, especially 
for narrow wells, except at low temperatures and high densities. Although a potential well of A = 1.5 is appropriate 
for many simple fluids, SW fluids with lower values of A may be of interest because the properties of certain colloidal 
suspensions are well reproduced by considering SW interactions with narrow potential widths. Therefore, as several 
theories for SW fluids have achieved a considerable accuracy for the equation of state and the pair correlation function 
of SW fluids, the constant volume excess heat capacity seems to be a suitable thermodynamic property to discriminate 
between them. In this context, we expect that our simulation data can stimulate other studies on the heat capacity 
of SW fluids of variable width and can be used to check the reliability of other approximations. 
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TABLE I: MC simulation data for Cy /Nk. The numbers enclosed between parentheses indicate the statistical uncertainty in 

the last decimal places. 





T* = 0.7 


T" = 1.0 


T* = 1.5 


I — z.U 


I — Z.O 








A = 1.1 






0.1 0.527(3) 


0.1731(8) 0.0580(2) 






0.2 


0.94(1) 


0.319(1) 


0.1109(3) 


0.0556(2) 




0.3 


1.25(2) 


0.431(3) 


0.1588(8) 






0.4 


1.40(2) 


0.525(5) 


0.1979(9) 


0.1033(4) 




0.5 


1.51(3) 


0.603(7) 


0.230(1) 






0.6 


1.50(3) 


0.630(7) 


0.252(2) 


0.1361(7) 




0.7 


1.42(3) 


0.612(7) 


0.261(3) 






0.8 


1.38(2) 


0.595(6) 


0.256(3) 


U.i42(ij 




0.9 


1.15(3) 


0.534(7) 


0.241(4) 












A = 1.2 






0.1 




0.367(2) 


0.1155(3) 


0.0559(1) 




0.2 


3.11(17) 


0.621(7) 


0.1988(6) 


0.0979(3) 




0.3 




0.77(1) 


0.2568(9) 


0.1288(6) 




0.4 


3.68(22) 


0.83(1) 


0.282(2) 


0.1463(6) 




0.5 




0.82(1) 


0.295(2) 


(j.ioZoyu ) 




0.6 


3.01(14) 


0.743(9) 


0.282(3) 


0.1459(8) 




0.7 




0.657(8) 


0.253(2) 


0.1351(7) 




0.8 


1.35(5) 


0.530(8) 


0.222(2) 


0.1zU9(9j 




0.9 




0.463(6) 


0.195(2) 


0.111(1) 










A = 1.3 






0.1 






0.1846(6) 


0.0864(2) 


0.0504(1) 


0.2 




1.17(1) 


0.303(2) 


0.1423(7) 


0.0834(3) 


0.3 






0.359(4) 


0.1720(5) 


0.1016(4) 


0.4 




1.35(2) 


0.363(2) 


0.178(1) 


0.1060(6) 


0.5 






0.340(3) 






0.6 




0.90(2) 


0.289(3) 


n 1 c: 1 1 \ 

(J.ioi^i j 


A AA1 O ^ 0\ 

\J.\J\)io[o ) 


0.7 






0.251(2) 


A 1'"»'~A/A\ 

0.13o0(9) 


0.0838(5) 


0.8 




0.512(8) 


0.226(2) 


0.iZD(i ) 


C "1 A O A "7 / \ 

U.UoU/ ) 


0.9 






0.219(2) 


0.1zz(lj 


A ATO / 1 \ 

0.078(1) 








A = 1.5 






0.1 






0.426(3) 


0.1724(4) 


0.0952(3) 


0.2 






0.705(5) 


0.263(2) 


0.1426(8) 


0.3 






0.719(9) 


0.277(2) 


0.1495(7) 


0.4 






0.563(5) 


0.239(2) 


0.1339(8) 


0.5 






0.401(6) 


0.190(2) 


0.1142(6) 


0.6 






0.295(2) 


0.161(1) 


0.1027(5) 


0.7 






0.270(2) 


0.1513(6) 


0.0977(6) 


0.8 






0.235(3) 


0.136(1) 


0.0894(9) 


0.9 






0.188(2) 


0.109(2) 


0.0716(6) 
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List of figures 




FIG. 1: Constant-volume excess heat capacity for a SW fluid witli A = 1.1 as a function of tiie reduced density p*. Circles: 
simulation data from Table^Jfor T* = 0.7, T* = 1.0, and T* = 1.5, respectively, from top to down. Squares: values obtained 
from Eq. (0J using the simulation data of F2 reported in j^fl]. Continuous curve: YAS model. Dashed curve: BH perturbation 
theory in the local compressibility approximation. Dotted curve: BH perturbation theory in the macroscopic compressibility 
approximation. 
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FIG. 2: As in Fig.^for A — 1.2, except that the temperatures are T* — 1.0, T* = 1.5, and T* = 2.0, respectively, from top to 
down. 
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FIG. 3: As in Fig. 0for A = 1.3, except that the temperatures are T* = 1.50, T* — 2.0, and T* = 2.5, respectively, from top 
to down. 




FIG. 4: As in Fig. 0for A = 1.5, except that the temperatures are T* = 1.50, T* — 2.0, and T* = 2.5, respectively, from top 
to down. 



